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Abstract 

A variational time discretization of anisotropic Willmore flow combined with a spatial dis¬ 
cretization via piecewise affine finite elements is presented. Here, both the energy and the metric 
underlying the gradient flow are anisotropic, which in particular ensures that Wulff shapes are 
invariant up to scaling under the gradient flow. In each time step of the gradient flow a nested 
optimization problem has to be solved. Thereby, an outer variational problem reflects the time 
discretization of the actual Willmore flow and involves an approximate anisotropic -distance 
between two consecutive time steps and a fully implicit approximation of the anisotropic Will¬ 
more energy. The anisotropic mean curvature needed to evaluate the energy integrand is replaced 
by the time discrete, approximate speed from an inner, fully implicit variational scheme for 
anisotropic mean curvature motion. To solve the nested optimization problem a Newton method 
for the associated Lagrangian is applied. Computational results for the evolution of curves un¬ 
derline the robustness of the new scheme, in particular with respect to large time steps. 


1 Introduction 

This paper generalizes a recently proposed variational time discretization [T| for isotropic Will¬ 
more flow to the corresponding anisotropic flow. Thereby, the anisotropic Willmore flow is de¬ 
fined as the gradient flow of the anisotropic Willmore energy with respect to the corresponding 
anisotropic L^-metric. 

The isotropic Willmore energy is given by w[x\ = ^ h^da, where x denotes the identity 
map and h the mean curvature on a surface M.. The isotropic L^-metric is given by (n, v)m = 

\v\‘^da , which is considered as a squared L^-distance of the surface M. being displaced with 
the vector field v from the non displaced surface M.. In the hypersurface case Willmore flow 
leads to a fourth order parabolic evolution problem, which defines for a given initial surface 
Ado a family of surfaces M.{t) for f > 0 with Ad(0) = Ado B9ll47l [30l . Applications of a 
minimization of the isotropic Willmore energy and the corresponding Willmore flow include 
the processing of edge sets in imaging l[36l [34l [5T] [T^l, geometry processing gSl |9l |8l |50|| 
and the mathematical treatment of biological membranes f2^ |46l |24l . Starting with work by 
Polden ll^ l4ll existence and regularity of Willmore flow was advanced in the last decade 

mMM- 

Now, in the context of Finsler geometry the classical area functional is replaced by the 
anisotropic area functional Si^[x\ = 'y{n)da with a local area weight 'y{n) depending on the 
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local surface orientation. Here, 7 is a positive, 1-homogeneous anisotropy function. In analogy 
to the isotropic case the anisotropic mean curvature h..y is defined as the L^-representation of 
the variation of the anisotropic area in the direction of normal variations of the surface and 
can be evaluated as h.y = div^^ (V 7 (n)). Hence, a possible first choice for an anisotropic 
Willmore functional is given by ^ da . Clarenz lITSll has shown that Wulff shapes are the 

only minimizers of this energy. Palmer studied variational problems involving anisotropic 
bending energies for surfaces with and without boundaries. Unfortunately, this energy definition 
does not imply the scale invariance property of Wulff shapes known for round spheres under 
isotropic Willmore flow. Indeed, any round sphere is a sfafionary point of the isotropic Willmore 
functional in M^. In a circle of radius Rq evolves under isotropic Willmore flow according 
fhe fhe ordinary differential equafion R = ^R~^- The counferparf of a round sphere in fhe 
anisofropic confexf is fhe Wulff shape as fhe unif ball wifh respect to the norm associated with 
the dual 7 * of the anisotropy 7 . But there is no such scaling law for the evolution of Wulff 
shapes under the above anisotropic variant of Willmore flow. 

To ensure full consisfency wifh fhe Finsler geomefry, one has fo adapf bofh fhe anisofropic 
energy and fhe anisofropic mefric as suggesfed in ll42ll (see Section [^. Indeed, we make use 
of fhe associated anisofropic mefric 7 * (v) (V 7 *) {v) ■ v'y{n)da (here only defined for 
v{x) 7 ^ 0 for all x G A4, cf. Section!^ for the general case), acting on a motion field v of 
fhe surface Ai wifh normal n. Furfhermore, we will use fhe anisofropic area weighf fo de¬ 
fine fhe anisofropic Willmore energy, i.e. Wj{x) = | h^ 7 (n)da. Then, if acfually fumed 

ouf fhaf Wulff shapes in acfually evolve according fo fhe same evolution law for radial 
paramefer valid for fhe evolution of circles under the isotropic flow. Recently, Bellettini & Mug- 
nai [T] investigated the first variation of this functional in the smooth case. Concerning the 
proper time and space discretization, this consistent choice of the anisotropic Willmore energy 
and the anisotropic metric on surface variations perfectly fits to the framework of the natural 
variational time discretization of geometric gradient flows. 

The finife elemenf approximafion of Willmore flow was firsl invesfigafed by Rusu ll44l based 
on a mixed mefhod for fhe surface paramefrizafion x and fhe mean curvafure vecfor hre as 
independenf variables, see also |[T6l for fhe application fo surface resforafion. In a level 
sef formulalion of Willmore flow was proposed. In fhe case of graph surfaces Deckelnick and 
Dziuk IfTSi were able fo prove convergence of a relafed space discrefe and fime confinuous 
scheme. Deckelnick and Schieweck established convergence of a conforming finite element ap¬ 
proximation for axial symmetric surfaces 1201 . In the case of the elastic flow of curves an error 
analysis was given by Dziuk and Deckelnick in ifT^ . An alternative scheme, which in partic¬ 
ular ensures a better distribution of nodes on the evolving surface was presented by Barrett, 
Garcke and Niirnberg ||2l|4l. Using discrete geometry calculus Bobenko and Schroder ifTOl sug¬ 
gested a discrete Willmore flow of friangular surfaces. The fime discrefizafion of fhe second 
order, anisofropic mean curvafure flow has been considered by Dziuk already in Il27ll2^ and he 
gave convergence resulfs for curves. Diewald 11211 has exfended fhe discrefizafion approach for 
isofropic Willmore flow of Rusu l44l fo some anisofropic varianf, for which Droske Il22l and 
Nemifz ll35l invesfigafed a level sef discrefizafion. 

Mosf of fhe above discrefizafion mefhods are based on some semi-implicif time discretiza¬ 
tion, which requires the solution of linear systems of equations at each time step. Thereby, the 
involved geometric differential operators are assembled on the surface from the previous time 
step. In the application one observes strong restrictions on the time step size. This shortcom¬ 
ing motivated the development of a new approach for the time discretization of Willmore flow 
in lU based on the following general concept for a variational time discretization of gradient 
flows: The gradienf flow on a (in general infinite dimensional) manifold with respect to an en- 
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ergy e[-] and a metric g on the manifold, is defined as the evolution problem x = —gradge[x] 
with initial data x^, where gradge[x] is the representation of the variation e'[x\ in the metric g, 
i.e. g{gTSLdge[x\X) = e'[x](C) for all infinitesimal variations Q of x. Now, one defines a time 
discrete family {x^)k=o,-- with the desired property x^ x{kT) for the given time step size r. 
To this end, one successively solves a sequence of variational problems, i.e. in time step k 

= arg min^ dist(x^, x)^ + 2 t e[x] , 

where dist(x^,x) = inf ^/gj(s)( 7 (s},j(s))ds denotes the Riemannian distance of 

,x\ V 

X from x^ on the manifold and r[x^,x] is the set of smooth curves 7 with 7(0) = x^ and 
7(1) = X. The striking observation for this abstract scheme is that one immediately obtains 
an energy estimate, i.e. e[x^'''^] + 27 dist(x*^, < e[x^]. In the context of geometric 

flows, this approach was studied by Luckhaus and Sturzenhecker ll^ leading to a fully implicit 
variational time discretization for mean curvature motion m BV and by Chambolle ifTTll . who 
reformulated this scheme in terms of a level set method and generalized it for the approximation 
of anisotropic mean curvature motion in [Ill . The time discretization for Willmore flow 
proposed in |T| builds upon this general paradigm. In this paper, we will show how to adapt the 
approach to the time discretization of the anisotropic Willmore flow which is fully consistent 
with Finsler geometry. 

The paper is organized as follows. In Section we briefly review the time discretization 
of isotropic Willmore flow. Building on these prerequisites the generalization to anisotropic 
Willmore flow is discussed in Section]^ Then, in Section|^we discuss a fully discrete numerical 
scheme based on piecewise affine finite elements on simplicial surface meshes. In Sectionj^the 
Lagrangian calculus from PDE constraint optimization is used to develop a suitable algorithm 
for the solution of the nested optimization problem to be solved in each time step. Finally, in 
Section [^computational results are presented. An appendix collects essential ingredients of the 
corresponding algorithm. 


2 Review of the time discretization of isotropic Will¬ 
more flow 

In this section we will briefly recall the nested time discretization of isotropic Willmore from 
lUl . We denote a hypersurface in by Af = Ai[y\. Here, y indicates a parametrization of 
M. and can also be considered as the identity map on Af parametrizing Af over itself. Then, 
the abstract variational time discretization of isotropic Willmore flow reads as follows: 

For given surface At [x^] with parametrization x^ and a time step r find a mapping x = x[x^] 
such that dist(Af [x^], Al[x])^ + r h^da —> min, where 

dist(Al[^], Ad[x])^ = — z)'^da is the squared -distance of surfaces M.[v] from the 

surface M.[z\,h. = h[x] is the mean curvature of Al[x], and da denotes the surface area 
of Af [x]. Now, we take into account that the mean curvature h = h[x] is the -gradient of the 
area functional on a surface Ad [x] and that mean curvature motion is the corresponding gradient 
flow. Thus, the mean curvature vector h[x]n[x] with n = n[x] denoting the normal on Ad[x] can 
be approximated by the discrete time derivative where y[x] is a suitable approximation 

of a single time step of the evolution of mean curvature motion with initial data x and time step 
size f. This time step itself can again be approximated using an (inner) variational scheme, i.e. 
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we define y[x] to be the minimizer of 


e,n[x,y] 



{y - xf + f\V M[x]y\'^da . 


(2.1) 


In fact, the corresponding Euler Lagrange equation is identical to the defining equation of the 
semi-implicit scheme for mean curvature motion proposed by Dziuk ||26| : 


0 = 



x)9 + fVM[x]y ■ '^M[x]dda 


( 2 . 2 ) 


Now, given y[x] as the minimizer of (^i for small f the functional | (vN^a;) 

approximation of the Willmore functional on Af [x]. This approximation is then used to define 
a variational scheme for a time step of the actual Willmore flow. To this end, we consider for 
given surface parametrization the functional 



x,y] 



x^yda + ^ {y — xfda , 
X J M.[x\ 


where we suppose y = y[x\ to be the minimizer of (^). To summarize, we obtain the following 
scheme for the kth time step of Willmore flow: 


Given an initial surface At[x’^] with parametrization x^ we define a sequence of surfaces 
At [x^] with parametrizations x^ for k = 1,... via the solution of the following sequence of 
nested variational problems 


= argmin^, e„„,[x^, X, y[x]], iv/zere (2.3) 

y[x\ = argmmj^e,„[x,?/] . (2.4) 


The inner variational problem p.4|) is quadratic, thus the resulting Euler-Lagrange equation 


( |2.2| ) is linear and we end up with a PDE constrained optimization problem to be solved in each 
time step. Eor more details we refer to |[T1. 


3 Nested time discretization for anisotropic Willmore 
flow 

Now, let us investigate the time discretization of anisotropic Willmore flow in the co-dimension 
one case. Here, we will in particular focus on the proper choice of energy and metric. We 
assume that 7 : [0, 00) is a positive, 1 -homogeneous (i.e. 7 (Ap) = |A| 7 (p) for all 

A G M,p G and sufficiently regular function, that satisfies the ellipticity condition 

l"{p)qq>Co\\qf \/p, q £ ,\\p\\ = l,p ■ q = 0 (3.1) 

for some positive constant cq and the Euclidean norm || • ||. As already mentioned 7 (n) represents 
the anisotropic area weight for a surface normal n. The isotropic case is recovered by choosing 
7 (-) = II • II. We define the dual function of 7 as 

7 *(x) := sup{(x,y)) | f G B-y} Vx G , 
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where denotes the unit Ball in the 7 -norm. The ellipticity assumption ensures that 7 ) 

and its dual space (M'^+^, 7 *) are uniformly convex Banach spaces and the duality map T : 
(M‘^+\ 7 *) ^ (M'^+\ 7 ), with 


nx) = ). 

is an odd single-valued bijective continuous map. More precisely T(0) = t),T{x) = {x)'V {x) 

for X 7 ^ 0, and = 7(C)^7(0 foi" C / 0- For details we refer to ll^ . The unit ball 

F := {x £ : 7 (x) < 1} in 7 ) is denoted the Frank diagram, the associated dual 

unit ball W := {x G : 7 *(x) < 1} is the corresponding Wulff shape. Wulff shapes are 

known to be solutions to the isoperimetric problem, that is dW minimizes the anisotropic area 
functional 

a.y[x] = / 7 (n[x])fia (3.2) 

J M. [x] 

(with 'y{n[x])da denoting the anisotropic area element) in the class of surfaces enclosing the 
same volume (cf. ifTTIl and the references therein). Now, based on the anisotropy 7 and its dual 
7 * we define an anisotropic distance dist-^ of a manifold Ai [y] from a manifold Af [x] by 


dist-y(Al[x], Ad[y])^ := f 'j*(y — x)‘^'y{n[x])da (3.3) 

J A4[x] 

for sufficiently regular x and y. The choice of the norm 7 * together with the anisotropic area 
weight 7 (n[x]) in ( |3.3| ) reflects the fact that the anisotropic area of the boundary of a convex 
body K C can be interpreted as 


a.y(dK) = lim 
' e-s>0 


\K + eW\ - \K\ 


where | • | denotes the usual Lebesgue volume in In particular, the underlying metric 

structure is dictated by the Wulff shape and its norm 7 * (see Q, l[42l and references therein). 

Based on these considerations let us first consider anisotropic mean curvature motion, which 
is defined as the gradient flow of the anisotropic surface area with respect to the above anisotropic 
metric. In this case the variational time discretization is associated with the minimization of 

dist-j.(Af[x], Af[y])^ + 2f f 'y{n[y])da (3.4) 

JMly] 

with respect to y for a given surface Af [x] and r > 0. Let us denote by y[x] the minimizer for 
given surface parameterization x. The Euler Lagrange equation for ( |3.4[ ) is given by 

0 = / T(y-x) • 6 l 7 (n[x])da + f(a^^[y], 6 ») 

JmIx] 

= f[ T ■ 9j{n[x])da + f{a^[y],e) (3.5) 

Jm[x] V ) 

for smooth test functions 0 : Af[x] —)■ Together with dty{kf ) ^ this reflects the 

weak formulation of anisotropic mean curvature motion given by 



T{dty) ■0-f{n[y])da = -{a^'[y],e) 


(3.6) 
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for a parametrization y and smooth test functions 9 defined on M. [y] (cf. Il42l i. Here, the varia¬ 
tion of the anisotropic area functional is given by 

{a^'[y],e)=l' h^[y]-p^ ■ e-f{n[y])da, 

JMly] l[n[y\) 

where h^[y] = (Prvj^y^{n^[y]) = div_yv([j^](V 7 (n[y])) denotes the anisotropic mean curvature 

with n^[y] = S/^{n[y]) (see HU). Thus, from ( |3.6[ ) we deduce that T{dty) = — or 

equivalently we achieve the strong formulation of anisotropic mean curvature motion 

dty = i^^[y\ := T-^ (^-h^[y] = -h^[y]V-f{n[y]). 

Indeed, as pointed out in H2ll the last equality holds due to the 1-homogeneity of 7 , i.e. 


7 



n[y] \ 
lin[y])) 


V 7 



-h^[y]V'y{n[y]). 


Next, we deal with the actual anisotropic Willmore flow and consider the anisotropic Will- 
more functional defined as follows for a parametrization x of Ad [x]: 

■= 7 / 7 *(K^[x])^ 7 (n[x])da. (3.7) 

^ JM[x] ^ JJvl[x] 

Here, we have used that the 1-homogeneity and V 7 (^) G dW for all ^ G imply 

= 7 * (-h7V7(n))^ = {-V'y{n)f = . 

Then the abstract variational time discretization of anisotropic Willmore flow reads as follows: 
Given Ad [x^] and time step r find a mapping x = x[x^] such that x minimizes 


dist-),(Ad[x^], Af [x])^ + r f 'y*{K^[x])‘^ 'y{n[x])da . (3.8) 

J M{x] 

As in the isotropic case, we will now replace the anisotropic mean curvature vector by the 
discrete speed extracted from a scheme for a single time step of anisotropic curvature flow ( |3.6| ). 
In explicit, 7 *( ^M ~^ )2 a suitable approximation of h^[x] = 7 *(«:-,,[x])^, where is the 

time discrete speed extracted from the variational time discretization of anisotropic curvature 
motion. Furthermore, we use the definition of the anisotropic distance measure in ( |3.3| ). Finally, 
based on this approximation we derive the actual time discretization of anisotropic Willmore 
flow. For a given surface parametrization x^ of the surface Ad [x*^] at a time step k we define the 
functionals 


e„„,[x^,x,y] := 

f 7*(x — x^)^ 7(n[x^])(ia-h 

^ / - 1 *{y - xf-i{n[x\)da 


J Ai[x^] 

'T Jm[x\ 

e.„[x,y] := 

/ 7*{y — x)‘^ 7(n[x])da -|- 2 f 

/ l{n[y])da, 


J M[x] 

' M[y\ 


and in analogy to the isotropic case above, we end up with the following fully nonlinear varia¬ 
tional time discretization of anisotropic Willmore flow: 
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Given an initial surface with parametrization we define a sequence of surfaces 

Ai with parametrizations for k = 1 ,... via the solution of the following sequence of 
nested variational problems 

x^~^^ = argmin^, a:, y[x]], where (3.9) 

y[x\ = arg miiiy e,„ [x, y] . (3.10) 


Different from the variational scheme for isotropic Willmore flow, the inner variational problem 
is no longer quadratic. It is worth to mention that this variational time discretization does not 
involve derivatives of the anisotropy. Nevertheless, as we will discuss below in the context of 
the actual computation, differentiation is required to run Newton methods for the associated 
Lagrangian functional. Indeed, for this we will need 7 , 7 * G C'^(M'^"'“^ \ {0}); moreover, unless 
( 7 *)^ G C^(M‘^+^) (which holds for 7 (p) = y/Ap ■ p with a symmetric positive definite matrix 
A), a regularization will be is required (see Section|^below). 

Let us conclude this section with a study of boundaries dW of two-dimensional Wulff 
shapes W moving under anisotropic Willmore flow in the plane. To this end consider the 
parametrization x : (0,T) x ^ M^, x{t, u) = of the boundary of a (rescaled) 

Wulff shape R{t)W. Using the results given in Il42l it is easily seen that x moves under anisotropic 
Willmore flow if R{t) solves the ODE 


R{t) 


1 

2R{tY ■ 


Hence, we observe that Wulff shapes expand in time like in the isotropic case (cf. [Til) with 
R{t) = ^ i?(0)^ + 2t . Next let us compare this with the time discrete evolution based on the 
proposed nested variational time discretization. We write x, y, x^ : 5*^ —)■ x{v) = RV'^{v), 

y{v) = RVj{i'), x^{v) = R^Vj{i'). Since 7 *(V 7 (z/)) = 1 we immediately derive 


eo„,[x^,x,y] = {R-R^)‘^ai^{x^) + ^{R-Rf‘aL^{x), 
ei„[x,y] = {R-Rfa^{x)+ 2fai^{y). 

Considering variations ye(i^) = (i? + in direction of the anisotropic normal we 

infer from the inner problem that 


(i? — R)ai^{x) + Xa.y(y) = 0. 

More precisely, since ay(y) = ■^a..y(x) due to the homogeneity property of 7 , we have that 

This, together with a.y(x) = -^aLy{x^), gives 

X, y] = a..^(x*^) (^{R - R’^f + , 

from which we deduce 

R-R^ _ 1 

T “ zrW 

Note that this is a slightly different time step scheme than the one reported for the isotropic case 
(7(') = II ■ ID in HI § 2 . 1 ]. This is due to fact that we use an implicit formulation of the inner 
problem as opposed to the linear equation p.l| ) in the scheme for isotropic Willmore flow (cf. 
Section 1^. 
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4 Finite element discretization in space 


Following the approach in |T| we now derive a suitable spatial discretization based on piecewise 
affine finite elements. This is in close correspondence to the surface finite element approach 
by Dziuk Il25l . To this end we consider simplicial meshes Af[X] as approximations of the 
hypersurfaces M.[x\ in i.e. polygonal curves for d = 1 and triangular surfaces for d = 2. 
Thereby, X is a parametrization of the simplicial mesh Ad[X] which is uniquely described 
by a vector X of vertex positions of the mesh. Here, and in what follows, we will always 
denote discrete quantities with upper case letters to distinguish them from the corresponding 
continuous quantities in lower case letters. Furthermore, a bar on top of a discrete function 
indicates the associated vector of nodal values, i.e. X = {Xi)i^j, where Xi = {X^, • • • , Xf~^^) 
is the coordinate vector of the ith vertex of the mesh and I denotes the index set of vertices. 
For d = 1 each element T is a line segment with nodes Xq and Xi (using local indices) and 
for d = 2 the elements T are planar triangles with vertices Xq, Xi, and X 2 and edge vectors 
Fq = X 2 — Xi, Fi = Xq — X 2 , and F 2 = Xi — Xq. Given a simplicial surface Ad[X], the 
associated piecewise affine finite element space is given by 

V[M[X]) := {U G C°(Af[X]) | U\t G Pi VT g M[X]} 


with the nodal basis denoted by Here, Vi is the space of affine functions on a simplex 

T. Thus, for U G V{Xi[X]) we obtain U = U = {U{Xi))i^i. Let us em¬ 

phasize, that the parametrization mapping X itself is considered as an element in V(Ad [X])'^"''^ 
and we recover the vector of nodes X = (Xj)ig/. 

With these algorithmic ingredients at hand we now can derive a fully discrete nested time 
discretization of anisotropic Willmore flow, as the spatially discrete counterpart of (3^i and 

drlol ): 

Given a discrete initial surface XI [X*^] with discrete parametrization X^ we compute a 
sequence of surfaces A4 [X^] with parametrizations X^ by solving the nested, finite dimensional 
variational problems 

X^+^ = argminjfgv(_v,[jYfc])d+i <f»,[X^,X,y[X]], where, (4.1) 

y[X] = argminygv(;^[jf])d+i £:,„[X,y]. (4.2) 


Here, the functionals and are straightforward spatially discrete counterpart of the func¬ 
tionals ei„[x, y] and y\ and defined by 

f,„[X,y] := [ 1 {j*{Y-Xf) j{N[X])da +2f [ 7 (X[y])da, 

Jm[x] Jm[y] 


^„„,[X^X,y] := [ i(7*(X-X^)2) 7 (X[X^])da 

[ I( 7 *(y-X) 2 ) y{N[X])da, 

Jm[x] 

where the nodal interpolation operator I renders the resulting scheme fully practical. To simplify 
the exposition, we introduce the discrete quadratic form X] = I ( 7 *(y)^) 7 (X[X])da 

(a nonlinear counter part of the quadratic form induced by the lumped mass matrix) and the dis¬ 
crete anisotropic area functional A-j,[y] = 7 (X[y])da, both of which are assembled 

from local contributions on simplices of the underlying simplicial grid Th- 


M^[Z, X] 


E 

T£Th 


id+l)\ 


Y, l*{ZT,f\y{RT[X]) 






(4.3) 
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= E hiRTlX]) ( 4 . 4 ) 

TeTh 

Here, Rt[X] = - Xt,o) for d = 1 and Rt[X] = {Xt,i - Xt,o) A (Xr ,2 - Xt,o) 

for d = 2. Hence, we can rewrite 

= M^[X-X’^,X^] + ^M^[Y-X,X], 

= M^[Y -X,X] + 2fA^[Y]. 

The necessary condition for Y[X] to be a minimizer of E^^[X, •] is given by the corresponding 
discrete Euler Lagrange equation 

0 = 5r£:i„[A,y[A]](0) = dzm^[Y - X,X](0) + 2fdYA^[Y]{Q) 
for all 0 e V{M[X]f+^. 

5 Optimization algorithm for the time steps 

In this section, the actual optimization algorithm for the nested, fully discrete variational prob¬ 
lem derived in Section|^is presented. Thereby, we apply a step size controlled Newton method 
(cf. P31 section 7) for the corresponding Lagrangian (cf. Nocedal & Wright OTIH . In our context 
the Lagrangian function for problem ( | 4 . 1 [ ), ( | 4 . 2 [ ) is given by 

C[X,Y,P]=£^^\X\X,Y]-dY£dX,Y]{P) 

for independent unknowns X,Y ^ and the Lagrange multiplier P g K(rf-ti)|/| (with a 

slight misuse of notation, we consider these unknowns as finite element function in the spaces 
V{Ai[X^]Y~^'^ and V(A4[X])‘^+^, respectively, or as the associated nodal vector in 
Lor an extensive discussion of the Lagrangian ansatz we refer to Il38ll . Now, we ask for crit¬ 
ical points {X,Y,P) of L. Indeed, 0 = 5p£[X,y,P](0) = y](0) is the Eu¬ 

ler Lagrange equation of the inner minimization problem with respect to Y for given X and 
0 = dyCiX ,Y, P]{Q) = X, y](0) — y](P, 0) is the defining equafion 

for fhe dual solufion P given Y as fhe solufion of fhe above Euler Lagrange equation. Linally, 
fhe Euler Lagrange equafion for fhe acfual consfrainf opfimizafion problem coincides wifh 

0 = dxC[X,Y,P]{e) = dxEoAX^XXK®) - ax9yfi„[X,y](P,0). 

Lor fhe gradienf of fhe Lagrangian £ we obfain 

/ dxSout — dxdySiniP) \ 
grad£= dySout - d^£i,{P) 

V -dyS^n J 

wifh 

5x£„„.[X^X,y](0) = dzM^[X-x^x’^m 

+^idxM-,[Y - X, X](0) - dzM^[Y - X, X](0)), 

9y£„„.[X^X,y](0) = ^dzM^lY -X,X]{e), 
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5x5yf.„[X,y](P,0) = -dlM^[Y-X,X]iP,e) + dxdzM^[Y-X,X]{P,e), 
d^£,^[X,Y]{P,Q) = dlM^[Y-X,X]{P,e) + 2fd^A^[Y]{P,Q). 


The Hessian of C, which is required to implement a Newton scheme, is given (in abbreviated 
form) by 

/ - dj^dyS^niP) dxdySout-dxd^£,„{P) -dxdy£, \ 

Hess£= dxdy£a,t - dxd^£n{P) d^£a,t - d^£in{P) 

V -dxdy£, -d\£, 0 / 


The different terms in Hess C are evaluated as follows: 


d‘x£aM{Q, 'h) 


dydx£oM{Q, 

<9y‘^^out(0) 

d\dy£,,{Q,^,'S) 

5x9y<?i„(0,^,H) 

af-£:.„( 0 ,T',H) 


= dlm^[X - X\ X"](0, M/) + ^ {djcM^lY - X, X](0, 
-2dzdxM^[Y - X,X](0, + 4M^[y - X,y](0, ^)) , 

= ^ (dzdxM .,[y - y, y] ( 0 , T-) - 9|M^[y - X, X]( 0 , ^)) , 

= ^dlM^[Y-X,X]iQ,^), 

= [y - X, X] (0 , t-, h) - dxdlM^ [y - x, x] ( 0 , h) 

-dxdlm^[Y - X, X](0, H, + dj,dzM^[Y - X, X] (0, T-, E), 

= -4M^[y - X, X] (0, ^,E)+ dxdlM^[Y - X, X] (0, ^,E), 

= dlM^[Y - X, X](0, M/, H) + 2f9f-A^[y](0, H). 


In the implementation of the proposed scheme it is convenient to directly treat the squared, dual 
anisotropy 7 *’^(.) := ( 7 *(-))^ the calculation of derivatives of the anisotropic functionals, 
which is particularly advantageous for anisotropies of the type 'y{p) = Ylk=i Vp ' GkP where 
the Gk are symmetric and positive definite (cf. Garcke et al. ||3l). The different terms of the 
gradient grad £ and the Hessian Hess C are in the usual way assembled from local contribution 
on simplices of the polygonal mesh. The required formulas are given in the Appendix. 


6 Numerical results 


In this section, we show applications of the proposed algorithm to the evolution of curves in 
under anisotropic Willmore flow. Beside anisotropies with ellipsoidal Wulff shapes we study 
regularized crystalline anisotropies 7 (-) = || • H^i and 7 (-) = || • based on a suitable 
regularization. A particular emphasis is on the verification of the robustness and stability of the 
proposed approach in particular for large time steps. Furthermore, we experimentally verify that 
Wulff shapes grow self-similar in time under the corresponding anisotropic Willmore flow. 

At first, we study anisotropies of the type 



2„2 


a^z 


1^1 


2„2 


+ ap2 


for given oi, 02 > 0. In that case the squared dual anisotropy function is given by 







A Nested Variational Time Discretization for Parametric Anisotropic Willmore Flow 


11 




Figure 1: The evolution of an unit eirele under isotropie Willmore flow is plotted on the left. For the 
eomputation we used as initial grid size h = 0.0981 resulting from 64 vertiees. Furthermore, t = h, 
f = and the resulting diserete eurves are shown for f = 0, lOr, 50r, lOOr, 500r. In the middle we 
display the evolution of an ellipse (with half axes 6 and 1) under anisotropie Willmore flow with 256 
elements and h = 0.0984. Here, we eonsider t = h, f = and display the approximate solutions 
for f = 0, lOr, 50r, lOOr, 500r. Next, the assoeiated L^-errors are plotted over time on the right, 
where the lower error eurve eorresponds to the evolution results on the left. 

Figure [T] compares the evolution of a circle of radius iio = 1 under isotropic Willmore flow for 
ai = 02 = 1 with the evolution of an ellipse with half axes oi = 6 and 02 = 1 under the corre¬ 
sponding anisotropic flow. As discussed in Section]^ in both cases the initial curve Ado expand 
in a self-similar fashion, i.e. Ad[a:(t)] = R{t)Aio with R{t) = i?Q -|- 2t for tq > 0. In Figure 
[^we plot the evolution of the error err{h) := \\Thx{t) — Xh{f)\\L 2 in time. Thereby, the L^-error 
is evaluated on the polygonal curve Xh{t) and Xh denotes the nodal interpolation of x{t) at the 
projected positions of the nodes of Xh{t) in direction Vj{n[xh{t)]). In Table[T]and|^we provide 
results on the experimental order of convergence eoc := log(err(/ii)/err(/i 2 ))/ log(/ii//i 2 ) for 
varying grid and time step size in case of the evolution of the circle and the ellipse. 

Now, we want to study crystalline anisotropies '}{■) = || • H^i and 7 (-) = || • ||^^. As 
already pointed out, even though the formulation of the scheme itself doesn’t explicitly need 
assumptions on the smoothness of 7 , the application of the optimization algorithm requires the 
computation of derivatives of 7 up to order 3. In fact, we use the following regularization: For a 
small parameter e > 0 we regularize the ^^-norm by 

1=1 

Since in the f“-norm equals a rotated and scaled £^-norm we use as regularization of the 
^ 00 -norm 

^ -b {zi + Z 2 )‘^ , - ^2)2 

V)- 2 + 2 

Figure 1^ shows the evolution of a sphere with respect to the regularized f“-norm under the 
associated anisotropy Willmore flow with anisotropy 7 (-) = || • H^i for s = 0.0001. Results 
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n 

h{t) 

L^-error 

(r = 7 = hi) 

eoc 

h{t) 

L^-error 
(r = r = ho) 

eoc 

4 

4.166e-l 

4.830e-3 


4.482e-l 

1.916e-2 


5 

2.096e-l 

1.328e-3 

1.879 

2.258e-l 

1.087e-2 

0.826 

6 

1.049e-l 

3.403e-4 

1.969 

1.132e-l 

5.804e-3 

0.909 

1 

5.249e-2 

8.561e-5 

1.992 

5.668e-2 

3.000e-3 

0.954 

8 

2.625e-2 

2.144e-5 

1.998 

2.836e-2 

1.525e-3 

0.977 


Table 1: The L^-error between the exaet solution of the self-similar evolution of eireles under Will- 
more flow and the diserete solution of the fully implieit variational time diseretization is plotted at 
time t = 0.1542 for a grid size h{t) (left) and t = 0.3927 (right). On the left we eonsider time step 
sizes r and f of the order of the squared spatial grid size ho at the initial time 0, whereas on the right 
both time step sizes are taken equal to the grid size. In both oases we have oonsidered 2"^ vertioes 
for the polygon, resulting in an initial grid size ho = . 


n 

h{t) 

L^-error 
(r = r = hi) 

eoc 

h{t) 

L^-error 

(r = f = ho) 

eoc 

5 

1.435e-r0 

1.648e-l 


1.214e+Q 

1.942e-l 


6 

6.487e-l 

3.476e-2 

1.960 

5.875e-l 

7.089e-2 

1.303 

7 

3.069e-l 

8.762e-3 

1.841 

2.842e-l 

3.424e-2 

1.002 

8 

1.525e-l 

2.182e-3 

1.987 

1.396e-l 

1.724e-3 

0.966 


Table 2: As in Table experimental orders of eonvergenee are reported, now for the self-similar 
evolution of the ellipses (with half axis 6 and 1) under anisotropio Willmore flow. Here, again 
polygons with 2” vertiees are oonsidered, equi-distributed along the initial ellipse with an initial 
grid size ho = On the left the error is evaluated at time t = 0.596576 and on the right at time 
t = 0.77238. 


on the self similar evolution of spheres with respect to the regularized f^-norm are depicted 
in Figure In these simulations, we use the analog regularization for the dual anisotropy 7 * 
required in the algorithm. 

Next, we generalize Willmore flow and replace the Willmore energy by the modified energy 

e^[x] := [ (l-hi + x] j{n[x])da, ( 6 . 1 ) 

JmIx] J 

with a second term given by the anisotropic area weighted with a constant A > 0. The in¬ 
corporation of this generalized energy in our computational approach is straightforward. The 
generalized flow combines expansive forcing with respect to the anisotropic Willmore flow of 
curves with contractive forcing due to the anisotropic mean curvature motion associated to the 
anisotropic area functional. Thus, for the generalized model we expect convergence to a limit 
shape given by a scaled Wulff shape, where the scaling depends on the factor A. Fig.|^shows the 
impact of the factor A on the evolution, whereas in Fig. we compare the evolution of different 
initial shapes under the generalized anisotropic Willmore flow for different anisotropies. 
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Figure 2: Evolution of the unit sphere with respect to the regularized £°°-norm under anisotropic 
Willmore flow for the anisotropy || ■ H^i with e = 0.0001. For this computation we consider 200 
vertices leading to an initial grid size ho = 0.04. Furthermore, t = ho and f = h^ and the resulting 
discrete curves are shown for t = 0, lOr, 50r, lOOr, 200r. 


Appendix 

Here, we eolleet the eomputational ingredients to evaluate the Lagrangian, its gradient and Hes¬ 
sian based on a standard loeal assembly proeedure. In the following for veetors x G and 
funetions / we use the notation fij{x) = and in analogy for higher order derivatives. 

Furthermore, for matriees A we use fk,ij{A) = and again in analogy for higher order 

derivatives. In faet, we ean restriet ourselves to the loeal funetionals 


i=0,...,d 


Mt,^[Z,X] = I 7(i?[A]), At,-,[A] = ^^iR[X]), (6.2) 


where we denote by Z = (Zq, ..., Z^) and X = (Xq, ..., Xd) the eorresponding veetors 
of simplex nodes in with eoordinate representation Zj = {Zjr)r=i,...,d+i and Xj = 

{Xjr)r=i,...,d+i- Here, i? is a mapping from to representing the 90° rotated edge 

veetor for d = 1 and the eross produet of edge veetors for d = 2, respeetively. For d = 1 we 

_ f X 2 _ ^12 \ 

obtain for the first derivatives of i?[X] = ( I with respeet to the entries {ij) with 


obtain for the first derivatives of i?[X] = 
i = 0,..., d and j = 1,..., d -|- 1 

Roi[X] =( M , R^o2[X] = ( I 


Xn — Xoi 


^-Oi[A] = (^ _1 J > RMX] =[o ) ^ ^MX] = [i ) ^ ^ 

Beeause of the linearity of d? for d = 1 all higher derivatives vanish. For d = 2 we have 


R[X] = - Xo„)(X2„ - Xov) 


i=l,2,3 


where e^uv is the Levi-Civita symbol (e^uv = ±1 if {ru, u, v) is a even/odd permutation of 
(1,2, 3) and 0 else). Thus, for tu = 1, 2, 3 we have 


Rw,js\X\ — ^ ^ (^wuv 'ioj)^su(X2v Xov) T (d2j doj)dgi,(Xi^ Xou)) ) 
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Figure 3: Evolution of the unit sphere with respect to the regularized £^-norm under anisotropic 
Willmore flow for the anisotropy 7 (-) = || • ||^oo. The parameters are ho = 0.0078, £ = 0.001, 
T = f = hi and curves are plotted at times t = 0, lOr, 50r, lOOr, 500r, lOOOr on the left and 
ho = 0.0283, £ = 0.0001, r = ho, f = hi, t = 0, lOr, 50r, lOOr, 200r, 275r in the middle. On the 
right the associated -errors are plotted over time, where the lower error curve corresponds to the 
evolution results on the left. 


R. 


'Wjs It 


[^] 


3 

{{5ij — Soj){ 62 l — 5oi)5su^tv + {^2] — Soj){6il — ^Ol)^sv^tu) , 

u,v=l 


and all third derivatives Ru^irjsit[X] vanish. Here j, I G {0,1, 2} refer to the local node and 
s,t £ {1,2,3} to the spacial component. Next we derive expressions for the derivatives of 
Q and \AA\ under the assumption that 7 , 7 * are sufficiently smooth in \ {0} ( thus 
/O): 




Figure 4: The impact of the parameter A is shown for the evolution of a circle to an ellipse with 
aspect ratio 4 : 1 (i.e. ai = 4 and 02 = 1). We evolve polygons with 160 vertices approximating the 
unit sphere as initial curve, ho = 0.0393 and r = f = 0.01, h = 0.000393. On the left A = 0.025 
and on the right A = 4. 
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Figure 5: The evolution of different initial shapes for different anisotropies is displayed. For all 
computations we use 100 vertices and choose A = 0.25. On the left we start with an ellipse with 
aspect ratio 4 : 1 under an isotropic flow with 7 (-) = || ■ || (ho = 0.1739, t = ho, f = h^) results 
are shows at t = 0,0.1739,0.5218,1.739,3.478,6.956,173.9. In the middle and on the right an 
ellipsoidal anisotropy with aspect ratio 2 : 1 is used (i.e. ai = 2, 02 = 1) in the first case (middle), we 
take as initial shape the unit sphere for the /^-norm (ho = 0.0566, r = f = 0.001 ho) and results are 
displayed at f = 0, 0.00017, 0.00085, 0.00169,0.006,0.056,0.251). In the second example (right), 
the initial shape is the unit sphere for the /°°-norm (ho = 0.08 and r = f = 0.01 ho) and results are 
depicted for t = 0, 0.0024, 0.008, 0.04, 0.08, 0.8,4.8. 


Derivatives of ^ 


dz,MT,^[Z,X] 
dz,,dz„MT,^[Z,X] 
dzit dzjs dzir Mt,7 

dx,MTM^,X] 


^^^jf(Z.)j(ElX]), 


dx,AMTM^,x] 

dx,,dz,MTM^,x] 
dxit dzjs dzir Mt,7 l^,N] 
dxit dxjs dzir Mt,7 l^,N] 


1 


(d+1)! 




> a=0,...,d 


f^J,t(ElXj)Et,irJslXj+ 7MR[X])Rt,irmRu,js[X]] , 
= 1 

m 

^^^jf(z,)j2xt(R[x])Rtjsm , 

^i^^^jf^(z,)Y,7,u(R[x])Ru,it[x] , 

' ' 11=1 


(d + 1) 


xf(Zi) 
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Y,iAR[x])Rvjsit[x] + 


Y, %vu{R[X])R,js[X]Ru,it[X] 


. V=1 


v,u=l 


Derivatives of ^ 


dY„AT,y[X] 

^YjAy^XtAX] 

dYitdYjAYir-XTAX] 


^Y^AR[x])RsMx], 

5=1 


d\ 


YlAR[X])Rt,irjs[X]+ Y l,tuiR[X])Rt,ir[X]Ru,js[X] 


t=l 


t,u=l 


^ lit lit 

Y^ 'y,vu{R[X])Rv^irjs[X]Ru,lt[X] + 'Y,1,v{R[X])Rv,irjslt[X] 

u,v=l 1 ?=! 

m 

+ 'y,vuw{R[X])Rv^ir[X]RuJs[X]R.w,lt[X] 

u,v,w=l 

m m 

+ Y^ 'y,vu{R[X])Rv,ir lt[X]Ru,js[X] + Y^ 7,vuiR[X])Rv^ir[X]RuJslt 

u,v=l u,v=l 
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